String-like Clusters and Cooperative Motion in a Model Glass-Forming Liquid 
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A large-scale molecular dynamics simulation is performed on a glass-forming Lennard- Jones mixture 
to determine the nature of dynamical heterogeneities which arise in this model fragile liquid. We 
observe that the most mobile particles exhibit a cooperative motion in the form of string-like paths 
("strings") whose mean length and radius of gyration increase as the liquid is cooled. The length 
distribution of the strings is found to be similar to that expected for the equilibrium polymerization 
of linear polymer chains. 
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Many glass-forming liquids at temperatures above 
their glass transition exhibit unusual dynamical proper- 
ties associated with incipient glass formation . In par- 
ticular, relaxation measurements (NMR, light scattering, 
dielectric, mechanical) on glass-forming liquids typically 
give relaxation functions which are non-exponential and 
these observations have been interpreted in terms of "dy- 
namic heterogeneity" in the liquid structure l^jdj. Mea- 
surements of the translational and rotational diffusion 
of tracer molecules in "glassy liquids" have shown 
the breakdown of continuum hydrodynamics and these 
deviations have also been interpreted in terms of liquid 
heterogeneity However, the typical optical clarity 
of glassy liquids, as well as neutron scattering measure- 
ments over a wide wavelength range, do not provide evi- 
dence for local density fluctuations which could be read- 
ily associated with the heterogeneities inferred from dy- 
namic measurements It therefore remains an open 
question whether glassy liquids contain heterogeneities, 
and if so, what form they take. 

The determination of the structure of glass- forming liq- 
uids requires high resolution measurements in space and 
time. Molecular dynamics (MD) simulations of glassy 
liquids provide a powerful adjunct to experimental stud- 
ies, since the molecular structure and motions can be 
examined in great detail. However, this computational 
approach has its own limitations since large, well equili- 
brated samples must be followed for long enough times 
to explore the dynamics of glassy liquids 0. It is also 
not obvious how this microscopic information should be 
best analyzed to reveal the nature of heterogeneities in 
these liquids. 

Previous MD simulations on glassy liquids have pro- 
vided evidence for their heterogeneity. For example, dy- 
namic heterogeneity has been illustrated by averaging 
the properties of glass-forming liquids over subregions to 
find that differences in these averages persisted for signif- 



icant times Is) . Simulation has also shown that the tracer 
diffusion of particles in glassy liquids deviates from the 
Stokes-Einstein relation at low temperature which 
has been interpreted in previous measurements as im- 
plying heterogeneity Q. Other MD studies have sought 
to examine heterogeneity by tracking the motion of in- 
dividual particles and have revealed ample evidence for 
dynamic heterogeneity in the sense that large dispersions 
in the particle mobilities in space and time were apparent 
JlOt . However, no quantitative description of the struc- 
tural nature of heterogeneity in glass-forming liquids has 
emerged from these studies. Here we return to this prob- 
lem. 

Extensive MD simulations of a glass-forming liquid in 
which the molecules interact through a Lennard- Jones 
(LJ) potential were performed to determine the existence 
and nature of heterogeneities in this model liquid [ pd] |. 
The system studied in Ref. pd| ] was a three-dimensional 
binary mixture of 8000 particles having relative concen- 
tration of 80% of type A and 20% of type B particles, 
where the size of the A particles was 10% larger than the 
B particles. The well depth e and the potential range a 
of the Lennard- Jones (L-J) potential. 
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and the relative particle concentration were chosen to 
suppress crystallization and phase separation p^ . Pre- 
vious studies of this system |13| have shown that it ex- 
hibits the properties of a "fragile" jlj] glass-forming liq- 
uid, and that the relaxation of density fluctuations found 
for this model liquid compares well with the predictions 
of mode-coupling theory. Ref. [0 reported evidence for 
dynamical heterogeneities in this model glassy liquid, but 
their microscopic structure and dynamics were not inves- 
tigated in detail. In this Report, we investigate the sim- 
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ulated liquid described in Ref. in order to determine 
the morphology of the heterogeneities that arise in this 
glassy liquid ||l5| . 

We report all quantities in dimensionless units involv- 
ing the L-J parameters of the A particles The sys- 
tem was equilibrated at 10 different temperatures T in 
the range [0.451, 0.550]. For reference, the mode-coupling 
temperature Tc was estimated as Tc — 0.435 ||l3|. The 
density varied from 1.09 particles per unit volume at the 
highest T to 1.19 at the lowest T simulated. Configura- 
tion histories for up to 4 x 10^ MD steps following equili- 
bration were stored for each run. Adopting argon values 
for the Lennard-Jones parameters of the larger spheres 
||17| gives an observational time of 26 ns. 
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tides of a dense liquid are localized in the cage created 
by surrounding particles, and hence G's(r, t) is also well 
approximated by G°{r,t) in the limit of small t. For in- 
termediate t, simulations of glassy liquids generally show 
substantial deviations of Gs {r, t) from G° (r, t) (Fig. m . 




FIG. 1. Self part of the van Hove correlation function 
Gs(r,t) vs. r at time t* when the non-Gaussian parameter 
a is a maximum (inset) The dashed line shows the ex- 

pected curve if the distribution of particle displacements was 
the Gaussian function G°{r,t) with (r^) equal to the simu- 
lation value at t* . The mobile particles are defined as those 
particles for which the particle displacement r exceeds r*. 

Ref. ||ll[] proposed a criterion for identifying particles 
of enhanced mobility, and showed that the most mobile 
particles in the liquid were spatially correlated. In this 
approach, particle mobility is assigned by comparing the 
self part of the van Hove correlation function Gs (r, t) with 
the Gaussian function G°{r,t) 
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describing the probability distribution of a Brownian par- 
ticle with displacement r from the origin. At long times, 
Brownian motion is eventually established in a liquid, 
and hence Gs {t, t) is well approximated by G" (r, t) in 
the limit of large t. G°{r,t) also describes a Brown- 
ian particle strongly localized by a harmonic potential, 
where (r^) is independent of t . At short times the par- 



FIG. 2. A snapshot of a large cluster of nearest-neighbor 
mobile particles at T = 0.451. The lighter particles give the 
initial particle positions, while the darker ones (shown con- 
nected to them via line segments) designate the same particles 
at a later time t* . 

Deviation from the limiting Gaussian G°{r,t) is conven- 
tionally measured by the "non-Gaussian parameter" a(t) 
f7| 



a{t) 



5(r2(t))2 



(3) 



which vanishes for G°{r,t). Fig. |^ shows Gs{r,t) for 
the coldest simulation (T = 0.451) at the time t* |l|] 
where a{t) has a maximum (see inset). Previous work has 
shown that dynamical heterogeneity is strongly exhibited 
at this characteristic time |0 . The long tail of Gg [r, t* ) 
relative to Gs{r,t*) shows that there exist particles with 
enhanced mobility relative to ordinary Brownian motion 
|]ll| . This point was also emphasized by Thirumalai and 
Mountain |^ . Following Ref. 1 1 1 1 , we classify a particle as 
"mobile" if its displacement r at t* exceeds the distance 
r* where Gs{r,t) equals G°s{r,t*) (see Fig. 0). 
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FIG. 3. Distinct van Hove correlation function G^{r,t) at 
various t. Inset: Fraction of mobile particles that at t* 
are within a distance 0.6 from a position occupied by another 
mobile particle at t = 0. 



proportional to the equilibrium radial distribution func- 
tion g{r). Gd{r,t) is then a time-dependent counterpart 
to g{r), which provides information about the degree of 
correlation in particle motions. For Brownian particles, 
the peaks of Gd{r,t), which are well developed at short 
times, decay toward unity as particle positions decorre- 
late through Brownian motion . The depletion of par- 
ticle density at small t in Gd{r,t) for r < 1 graduaUy in- 
creases at later t as other particles move into the position 
previously occupied by the test particle. If the position 
of the original test particle is occupied with high proba- 
bility by another mobile particle at later times, (r, t) 
should develop a strong maximum about r = 0. This 
is shown in Fig. ||, which presents G™(r, t) at different 
t for the lowest T. At t > the height of the nn peak 
decreases, and a peak at r = develops and grows with 
time, reaching its maximum value near t ^ t* . Note that 
no significant peak exists for the Gd{r, t) between mobile 
and other particles. 



To analyze the geometrical configuration and dynami- 
cal motions of these mobile particles, we define nearest- 
neighbor (nn) mobile particles as those separated by a 
distance less than that of the first minimum in the radial 
distribution function. We can thereby identify connected 
clusters of nn mobile particles. A snapshot of a large clus- 
ter of nn mobile particles in the coldest run is shown in 
Fig. Ij. The size distribution for these clusters is broad, 
typically ranging from 1 to about 200 particles at the low- 
est T. Note that many of the mobile particles are within 
large clusters at low T, since the number of particles in 
the largest cluster is of the same order of magnitude as 
the total number of mobile particles. It should also be ap- 
preciated that this is a projection of a three-dimensional 
structure, and consequently its diffuse, string-like nature 
is somewhat obscured. The average number Nm of mo- 
bile nearest neighbors of a mobile particle increases from 
TV™ = 1.87 at the highest T to = 2.43 at the low- 
est T, consistent with the string-like appearance of the 
clusters (cf. Fig. ||. In comparison, the average coor- 
dination number of an atom in this system is approxi- 
mately 13. The spatial extent of the cluster in Fig. || is 
about 15 molecular units, which corresponds to ^ 3.3 nm 
for orthoterphenyl, for which the hydrodynamic radius is 
0.21 nm ig. 

An examination of the particle motions indicates a cor- 
related motion of particles along quasi-one-dimensional 
paths. To quantify the nature of this correlated motion, 
we calculate the distinct van Hove correlation function 
G™(r, i) for the mobile particles, which describes the 
time-dependent density of particles in the vicinity of an 
arbitrarily chosen test particle at t = 0. The term "dis- 
tinct" indicates that the test particle is excluded from the 
calculation of the density, and at short times Gd{r,t) is 
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FIG. 4. Contour plot of the probability distribution of the 
vector 5'ri{t*) in the plane defined by 5fi{t*) and Tij. The 
vector fij defines the positive x axis. The increment between 
any two contours is 5% of the total range. The distribution 
vanishes by definition in the circular central region. 

The inset of Fig. ^ shows the fraction (j) of mobile par- 
ticles that are within a distance 5 = 0.6 at t* from a po- 
sition occupied by another mobile particle at t = [p^ . 
The T-dependence of 0, as well as that of the height of 
the peak at r = 0, indicates that the tendency for a mo- 
bile particle to move toward a position previously occu- 
pied by another mobile particle increases with decreasing 
T . Collective motion was recognized in previous simula- 
tions of glass formation, but the effect was not studied 
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quantitatively. See e.g., |20|. 

The particular string-like character of the collective 
particle motion is demonstrated by examining the an- 
gular correlation between the motion of two neighbor- 
ing mobile particles. Fig. ^ shows a contour plot of 
the probability distribution of the displacement vector 
Sri{t*) — nit*) — ri{0), projected onto the plane defined 
by 5ri{t*) and = fi{0) — rj{0), where j denotes a 
mobile particle within the nn shell of a mobile particle 
i. The vector fij points in the direction of the positive x 
axis. The contours divide the range of values of the prob- 
ability into 20 intervals so that the probability inside the 
innermost contour centered on (0.7, 0.0) is at least 20 
times higher than the probability outside the outermost 
contour. The distribution has been calculated using all 
mobile particles that have at least one other mobile par- 
ticle in their nn shell. If more than one mobile nn particle 
exists in the proximity of a reference mobile particle, then 
the averaging for the data in Fig. ^ includes all distinct 
pairs of nn mobile particles. In the absence of any cor- 
relation between 6ri(t*) and uj, the contour plot should 
have a rotational symmetry about the origin. Indeed, 
a contour plot calculated in a similar way between other 
(not mobile) particles is rotationally symmetric. Instead, 
the pronounced asymmetry of the contour plot in Fig. ^ 
shows that it is much more probable for a mobile particle 
to move in the direction of another mobile particle than 
in any other direction. It is this dynamical correlation 
that appears to cause the string-like clustering of mobile 
particles. 
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FIG. 5. A large loop of 13 mobile particles exhibiting corre- 
lated string-like motion at T = 0.451. Line segments connect 
identical particles at successive times as in Fig. 2. 

An animation of the cluster of nn mobile particles 
shown in Fig. || reveals that not all of the mobile particles 
are involved in string-like collective motion. These ob- 
servations suggest that we separately examine nn mobile 



particles which move collectively in string-like paths. It is 
anticipated that these cooperatively rearranging string- 
like clusters could have a large influence on momentum 
pl| , p2t and mass transport |^ in glassy liquids, so we 
examine these structures in detail. The potential signifi- 
cance of cooperative motion in glass-forming liquids has 
long been appreciated [ p3| . 

We define "strings" by connecting two mobile parti- 
cles i and j if min[|f* (r) - rj{0)\, |f,(0) - rjit*)\] < 0.6. 
This condition implies that one of the mobile particles 
has moved, and a second mobile particle has occupied its 
position. Fig. ^ shows a representative string which 
has a notable loop-like form. Most of the strings have 
free ends, but loops occur with modest frequency. The 
strings have a smaller average size than the nn mobile 
particle clusters. For example, the cluster in Fig. || con- 
sists of about 10 substituent strings. The average length 
and size of the strings increase as the liquid approaches 
the glass transition and their mass n distribution P{n) 
is shown in Fig. ^ for three different T. The semi-log 
plot reveals that P{n) is approximately exponential, i.e. 
P{n) ~ exp(— n/(n)), where the average "string length" 
(n) increases with decreasing T (see inset). For small n, 
the average radius of gyration Rg of the strings appears 
to grow linearly with (n), while for large n, an exponent 
near 1/2 is observed, similar to discrete random walk 
chains. Observations similar to Fig. ^ have been found in 
studies of the equilibrium polymerization of linear chain 
polymers [p5[. 
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FIG. 6. Probability distribution P{n) of string lengths n 
for various T. Inset: First moment (n) of P{n) vs. T. 

By analogy with equilibrium polymerization [p5| , the 
fraction of mobile particles involved in strings defines an 
"order parameter" for describing the changing popula- 
tion of mobile particles. Since a particle exchange event 
involves at least two particles, we count the fraction of 
mobile particles involved in strings of length n > 2. This 



4 



fraction increases from 0.51 for the highest T studied to 
0.75 for the lowest T, although the total fraction of all 
particles that are mobile (0.055 ± 0.005) changes little 
over the T range investigated. That is, an increasing 
fraction of mobile particles participates in strings as the 
system is cooled. 

It should be noted that the occurrence of such dynam- 
ical heterogeneities does not conflict with recent neutron 
scattering observations Q which excluded the occurrence 
of compact clustering in propylene glycol, a model organic 
glass former. It is well known from previous neutron scat- 
tering studies on equilibrium polymerization in, e.g., liq- 
uid sulphur that detecting such string-like structures 
by neutron scattering is virtually impossible due to the 
small contrast between the chains and the unpolymerized 
species. This difficulty should exist for other polymeriz- 
ing liquids as well, e.g. Te and Se [p7| . 

We can make some comparisons with other well- 
studied physical problems where string-like structures 
have been predicted and observed. The tendency to form 
such structures, which involves a breaking of the local 
rotational symmetry, is known to occur in, e.g., liquid 
crystals and in many high energy physics contexts 
p9| , ^ . This type of spontaneous symmetry breaking has 
been suggested previously for glassy liquids and glasses 
pl| , but previous investigations have not emphasized the 
emergence of clusters which arise from correlated particle 
motion as in the present study p2[ . 

The observation of string-like heterogeneities in this 
cold, dense, equilibrium Lennard-Jones liquid naturally 
leads to many other questions about glassy liquids. For 
example, do strings exist in other model glass-formers? 
Can they be observed in real liquids? What is the na- 
ture of string-string interactions and how do these inter- 
actions relate to the occurrence of the glass transition? 
Can the characteristic properties of glasses at low T (spe- 
cific heat linear in T, slowly varying thermal conductivity, 
etc.) be explained by the present observation of string 
excitations? Are the excitations similar in "strong" and 
"fragile" liquids ||2^? Does the occurrence of physical ag- 
ing in glasses reflect a structural evolution of the string 
excitations? What effect does confinement have on the 
dynamics of the strings and how does this relate to ma- 
terial properties? The study of the geometrical form of 
dynamical heterogeneities in glassy liquids promises to 
lead to an improved understanding of many aspects of 
glass formation. 
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meinschaft under SFB 262. 



[1] See, e.g. the recent special issue of Science [297, 1945 
(1995)] devoted to problems and properties of glasses. 

[2] In the present work we use the term "glassy liquid" to 
refer to a liquid above the glass transition temperature 
Tg, but below the temperature where transient localiza- 
tion of the particle motion ("caging") first occurs. We 
avoid the term "supercooled" since it is not clear that 
this liquid can crystallize. Certainly, some liquids cannot 
crystallize (cf. E. A. DiMarzio, in Comp. Mat. Sci. 4, 317 

(1995) .) 

[3] K. Schmidt-Rohr and H.W. Spiess, Phys. Rev. Lett. 66, 
3020 (1991); J. Leisen, K. Schmidt-Rohr and H. W. 
Spiess, J. Non-Crystalline Solids 172-174 737 (1994); A. 
Heuer, M. Wilhelm, H. Zimmermann and H.W. Spiess, 
Phys. Rev. Lett. 95, 2851 (1995); C. T. Moynihan and 
J. Schroeder, J. Non-Crystalline Solids 160, 160 (1993). 

[4] F. Fujara, B. Ceil, H. Sillescu and G. Fleischer, Z. Physik 
B 88, 195 (1992). 

[5] M.T. Cicerone, F.R. Blackburn and M.D. Ediger, J. 
Chem. Phys. 102, 471 (1995); M.T. Cicerone and M.D. 
Ediger, J. Chem. Phys. 103, 5684 (1995); F. R. Black- 
burn, M. T. Cicerone, G. Hietpas, P. A. Wagner, M. D. 
Ediger, J. Non-Crystalline Solids 172-174, 256 (1994); 
L. Andreozzi, A. Di Schino, M. Giordano and D. Lep- 
orini, Europhysics Letters 29 (1997). See also [^. 

[6] R.L. Leheny, N. Menon, S. R. Nagel, D. L. Price, K. 
Suzuya and P. Thiyagarajan, J. Chem Phys. 105, 7783 

(1996) . 

[7] W. Kob, p. 1 in Vol. Ill of "Annual Reviews of Compu- 
tational Physics", Ed. D. StaufFer (World Scientific, Sin- 
gapore, 1995); S.J. Plimpton, Comp. Mat. Sci. 4, 361, 
(1995). 

[8] D. Thirumalai, R. D. Mountain and T. R. Kirkpatrick 

Phys. Rev. A 39, 3563 (1989). 
[9] D. Thirumalai and R. D. Mountain, Phys. Rev. E 47, 

479 (1993). 

[10] M. Hurley and P. Harrowell, Phys. Rev. E 52, 1694 
(1995); A. I. Mel'cuk, R. A. Ramos, H. Gould, W. Klein 
and R. Mountain, Phys. Rev. Lett. 75, 2522 (1995); T. 
Muranaka, Y. Hiwatari, Phys. Rev. E 51, R2735 (1995); 
R. D. Mountain, in "Supercooled Liquids: Advances and 
Applications", Proceedings of American Chemical Soci- 
ety, in press. 

W. Kob. C. Donati. S.J. Pl impton, P.H. Poole and S.C. 
Glotzer, |gond-mat/9706075 . 

The L-J parameters e^^ and ctc/j for a, /3 € {A, B} are 
given by: eaa = 1-0, gaa = 1.0, eas = 1.5, oab = 0.8, 
egg — 0.5, and gbb = 0.88. Note the presence of the sub- 
stantial A-B interparticle interaction which counterbal- 
ances the tendency for the size ratio oaaIcsbb to cause 
phase separation. 

W.Kob and H.C. Andersen, Phys. Rev. E 51, 4626 
(1995); W. Kob and S.J. Plimpton, unpublished. 
C.A. Angell, J. Non-Cryst. Solids 131-133, 13 (1991). 
The simulations discussed here were run on the massively 
parallel nCUBE 2 and Intel Paragon supercomputers at 
Sandia National Laboratories. In aggregate the equilibra- 
tion and statistics-gathering runs required approximately 
300 hours of CPU time on 256 Paragon nodes (i860 pro- 
cessors) and generated about 10 Gbytes of data. 
[16] The length is defined in units of oaa, temperature T in 



[11] 
[12] 



[13] 

[14] 
[15] 



5 



units of eAA/kB, where fcs is the Boltzmann's constant, 
and time t in units of ^ a\jjn/ eAA- The mass m of both 
types of particles are taken to be equaL 

[17] A. Rahman, Phys. Rev. A 136, 405 (1964); K. Skold, J. 
M. Rowe, G. Ostrowski and P. D. Randolph, Phys. Rev. 
A 6, 1107 (1972). 

[18] The time t* corresponds to the end of the transient 
plateau in {r^{t)) where the particle motion crosses over 
from a localized to asymptotic diffusive behavior. This 
time was shown to be related to the cluster lifetime in 
Ref. 0. 

[19] The qualitative results are insensitive to the choice of S, 
provided that S is chosen smaller than (jaa- 

[20] J.L. Barrat, J.N. Roux, J.P. Hansen, Chem. Phys. 
149, 197 (1990); H. Miyagawa, Y. Hiwatari, B. Bernu, 
J.P.Hansen, J. Chem. Phys. 86, 3879 (1988).; G. Wahn- 
strom, Phys. Rev. A 44, 3752 (1991); H. Teichler, 
preprint. 

[21] S. P. Chen, T. Egami and V. Vitek, Phys. Rev. 37, 2440 
(1988). 

[22] Douglas and Hubbard introduced a model of stress re- 
laxation in fragile glass-forming liquids based on the 
assumption of string excitations having an exponential 
length distribution, and found that many aspects of glass 
stress relaxation could be rationalized by their model. 
They also predicted that the structures dominating the 
viscoelasticity of "strong" glassformers should be sheet- 
like. J.F. Douglas and J.B. Hubbard, Macromolecules 24, 
3163 (1991); J.F. Douglas, Comp. Mat. Sci. 4, 292 (1995). 

[23] G. Adam and J.H. Gibbs, J. Chem. Phys. 43, 139 (1965); 
J. Jackie, A. Kronig, J.Phys.: Cond. Matter 6, 7633 
(1994). 

[24] For an animation of t he stringlike, correlated moti on of 
the loop in Fig. ^ see littp://www. ctcms.nist.gov/ ~do- 



nati/movie.html. 
[25] Y. Rouauh and A. Milchev, Phys. Rev. E 51, 5905 

(1995); S.C. Greer, Adv. Chem. Phys. 94, 261 (1996). 
[26] R. Bellissent, L. Descotes, P. Pfeuty, J. Phys.: Cond. 

Matter 6, A211 (1994). 
[27] J. Hafner, J. Phys. : Cond. Matter 2, 1271 (1990); 

G. Fairre and J.L. Gardissat, Macromolecules 19, 1988 

(1986). 

[28] I. Chuang, R. Durrer, N. Turok, and B. Yurke, Science 

251, 1336 (1991). 
[29] M. Stone, P.R. Thomas, Phys. Rev. Lett. 41, 351 (1978); 

T. Banks, R. Myerson and J. Kogut, Nucl. Phys. B 

129, 493 (1977); A.M. Polyakov, Nucl. Phys. B 120, 429 

(1997). 

[30] M.B. Hindmarsh and T.W. Kibble, Rep. Progr. Phys. 58, 
477 (1995). 

[31] N. Rivier and D.M. Duffy, J. Phys. (France) 43, 293 
(1982); Phil. Mag. A 40, 859 (1979); G. Venkataraman, 
Contemp. Phys. 26, 579 (1985); ibid. 27, 3 (1986). 

[32] String-like excitations involving correlated atomic mo- 
tions in liquid ^f/e were suggested long ago by Feynman 
(R.P. Feynman, Phys. Rev. 90, 1116 (1953); ibid. 1291 
(1953); ibid. 94, 262 (1954)) in his model of superfluids. A 
path-integral Monte Carlo implementation of Feynman's 
model can be found in D.M. Ceperley and E.L. Pollock, 
Phys. Rev. B 39, 2084 (1989). 



6 



